
#ifndef _EULER_H_
#define _EULER_H_

namespace Tuna {

template <typename T> class Euler;

template <typename T>
struct Typeinfo <Euler<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};


//template<typename Interpolator_t, int Dim>
template<typename Interpolator_t>
class Euler 
{
  typedef typename Typeinfo<Interpolator_t>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Interpolator_t>::Dim>::tiny floatTinyArray_t;
  typedef typename TunaArray<int, Typeinfo<Interpolator_t>::Dim>::tiny intTinyArray_t;
  typedef typename TunaArray<prec_t, Typeinfo<Interpolator_t>::Dim>::huge ScalarField_t;
///  Traits to define the type of arrays to be used.
///
///  floatTinyArray_t : Array of floats mostly used for deltas of meshes
///
/*
    typedef typename TunaArray<double, Dim>::tiny floatTinyArray_t;
    typedef typename TunaArray<int, Dim>::tiny intTinyArray_t;
    typedef typename TunaArray<double, Dim>::huge ScalarField_t;
*/
 public:
    Euler() { } ;
    ~Euler() { } ;
    
/*!
 * This function make the integration of the ODE \f$ \f$ to calculate
 * the new position of the particle. The result is stored in \c newPos.
 * \param oldPos Actual position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 * \param u X-component of the velocity.
 * \param v Y-component of the velocity.
 * \param dt Time step.
 */
    inline floatTinyArray_t integrate(floatTinyArray_t &oldPos,
				      const floatTinyArray_t &deltas,
				      ScalarField_t &u,
				      ScalarField_t &v,
				      prec_t dt)

	{
	  UVW = i_object.interpolate(oldPos, deltas, u, v); /// Interpolate de velocity
	  newPos(0) = oldPos(0) + dt * UVW(0); /// Use Euler to update the position.
	  newPos(1) = oldPos(1) + dt * UVW(1); /// in each coordinate.
	    return newPos;
	}

    inline floatTinyArray_t integrate(floatTinyArray_t &oldPos,
				      const floatTinyArray_t &deltas,
				      ScalarField_t &u,
				      ScalarField_t &v,
				      ScalarField_t &w,
				      prec_t dt)
	{
	    UVW = i_object.interpolate(oldPos, deltas, u, v, w);
	    newPos(0) = oldPos(0) + dt * UVW(0);
	    newPos(1) = oldPos(1) + dt * UVW(1);
	    newPos(2) = oldPos(2) + dt * UVW(2);
	    return newPos;
	}   
    
 private:
    floatTinyArray_t UVW;
    floatTinyArray_t newPos;
    Interpolator_t i_object;
};

} // Tuna namespace

#endif // _EULER_H_
